Areal kernel demonstration

Code to fit centroid kernel

rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())

sf <- mw

centroid_distance <- function(sf) {
  cent <- sf::st_centroid(sf)
  D <- sf::st_distance(cent, cent)
  return(D)
}

D <- centroid_distance(sf)

## Warning in st_centroid.sf(sf): st_centroid assumes attributes are constant over geometries of x

dat <- list(n = nrow(sf),
            y = sf$y,
            m = sf$n_obs,
            mu = rep(0, nrow(sf)),
            D = D)

nsim_warm <- 100
nsim_iter <- 1000

fit_centroid <- rstan::stan(
  "centroid.stan",
  data = dat,
  warmup = nsim_warm,
  iter = nsim_iter
)

## Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#bulk-ess

## Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#tail-ess

Code to fit integrated kernel

sf <- mw

# The number of integration points within each area
L <- 10

# The method for selecting integration points
type <- "hexagonal"

# The number of areas
n <- nrow(sf)

# Draw L samples from each area according to method type
samples <- sf::st_sample(sf, type = type, size = rep(L, n))

plot_samples <- function(samples){
  ggplot(sf) +
    geom_sf(fill = "lightgrey") +
    geom_sf(data = samples, alpha = 0.5, shape = 4) +
    labs(x = "Longitude", y = "Latitude") +
    theme_minimal() +
    labs(fill = "") +
    theme(axis.title = element_blank(),
          axis.text = element_blank(),
          axis.ticks = element_blank())
}

plot_samples(samples)

# Construct an (L * n) * (L * n) matrix containing the Euclidean distance between each sample
# Note that this ^ will not be exact for some settings of type (hexagonal, regular)
S <- sf::st_distance(samples, samples)
dim(S)

## [1] 273 273

# Data structure for unequal number of points in each area
sample_index <- sf::st_intersects(sf, samples)
sample_lengths <- lengths(sample_index)
start_index <- sapply(sample_index, function(x) x[1])

dat <- list(n = nrow(sf),
            y = sf$y,
            m = sf$n_obs,
            mu = rep(0, nrow(sf)),
            sample_lengths = sample_lengths,
            total_samples = sum(sample_lengths),
            start_index = start_index,
            S = S)

fit_integrated <- rstan::stan(
  "integrated.stan",
  data = dat,
  warmup = nsim_warm,
  iter = nsim_iter
)

## Warning: There were 1 divergent transitions after warmup. See
## https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## to find out why this is a problem and how to eliminate them.

## Warning: Examine the pairs() plot to diagnose sampling problems

## Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#bulk-ess

Comparison

Sampling speed

samples_per_second <- function(fit) {
  # Outside of warm-up phase
  times <- rstan::get_elapsed_time(fit)
  # Number of samples in first chain
  nsim_iter <- fit@stan_args[[1]]$iter
  nsim_iter / mean(times[, 2])
}
samples_per_second(fit_centroid)

## [1] 185.3997

samples_per_second(fit_integrated)

## [1] 3.562621

Convergence analysis

Trace-plots:

bayesplot::mcmc_trace(fit_centroid, pars = "l")

bayesplot::mcmc_trace(fit_centroid, pars = "beta_0")

bayesplot::mcmc_trace(fit_centroid, pars = "rho[28]") # For example

bayesplot::mcmc_trace(fit_integrated, pars = "l")

bayesplot::mcmc_trace(fit_integrated, pars = "beta_0")

bayesplot::mcmc_trace(fit_integrated, pars = "rho[28]")

Rank histograms:

bayesplot::mcmc_rank_overlay(fit_centroid, pars = "l")

bayesplot::mcmc_rank_overlay(fit_centroid, pars = "beta_0")

bayesplot::mcmc_rank_overlay(fit_centroid, pars = "rho[28]")

bayesplot::mcmc_rank_overlay(fit_integrated, pars = "l")

bayesplot::mcmc_rank_overlay(fit_integrated, pars = "beta_0")

bayesplot::mcmc_rank_overlay(fit_integrated, pars = "rho[28]")

All Rhat < 1.1 necessary but not sufficient:

which(rstan::summary(fit_centroid)[["summary"]][, "Rhat"] > 1.1)

## named integer(0)

which(rstan::summary(fit_integrated)[["summary"]][, "Rhat"] > 1.1)

## named integer(0)

Fitted values comparison:

fitted_centroid <- summary(fit_centroid)$summary[paste0("rho[", 1:28, "]"), "mean"]
fitted_integrated <- summary(fit_integrated)$summary[paste0("rho[", 1:28, "]"), "mean"]
plot(fitted_centroid, fitted_integrated)